1 SymPortal pre-submission

Submitting ITS2 samples to Symportal

Prior to submission: retaining only paired reads that match ITS2 primer (we pooled with 16S during sequencing)

# *note: most of this was written by Dr. Carly D. Kenkel

# in Terminal home directory:
# following instructions of installing BBtools from https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/installation-guide/
# 1. download BBMap package, sftp to installation directory
# 2. untar: 
tar -xvzf BBMap_(version).tar.gz
# 3. test package:
cd bbmap
~/bin/bbmap/stats.sh in=~/bin/bbmap/resources/phix174_ill.ref.fa.gz

# my adaptors, which I saved as "adaptors.fasta"
# >forward
# AATGATACGGCGACCAC
# >forwardrc
# GTGGTCGCCGTATCATT
# >reverse
# CAAGCAGAAGACGGCATAC
# >reverserc
# GTATGCCGTCTTCTGCTTG

# primers for ITS2:
# >forward
# GTGAATTGCAGAACTCCGTG
# >reverse
# CCTCCGCTTACTTATATGCTT

# making a sample list based on the first phrase before the underscore in the .fastq name
ls *R1_001.fastq | cut -d '_' -f 2 > samples.list

# cuts off the extra words in the .fastq files
for file in $(cat samples.list); do  mv ${file}_*R1*.fastq ${file}_R1.fastq; mv ${file}_*R2*.fastq ${file}_R2.fastq; done 

# gets rid of reads that still have the adaptor sequence
# [I stopped doing this 'cause it never happened]
#for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq ref=adaptors.fasta out1=${file}_R1_NoIll.fastq out2=${file}_R2_NoIll.fastq; done &>bbduk_NoIll.log

##getting rid of first 4 bases (degenerate primers created them)
for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq ftl=4 out1=${file}_R1_No4N.fastq out2=${file}_R2_No4N.fastq; done &>bbduk_No4N.log

# only keeping reads that start with the primer
for file in $(cat samples.list); do ~/bin/bbmap/bbduk.sh in1=${file}_R1.fastq in2=${file}_R2.fastq k=15 restrictleft=21 literal=GTGAATTGCAGAACTCCGTG,CCTCCGCTTACTTATATGCTT outm1=${file}_R1_ITS.fastq outu1=${file}_R1_check.fastq outm2=${file}_R2_ITS.fastq outu2=${file}_R2_check.fastq; done &>bbduk_ITS.log
# higher k = more reads removed, but can't surpass k=20 or 21

# renamed them to the shorter version again after checking
for file in $(cat samples.list); do  mv ${file}_*R1*.fastq ${file}_R1.fastq; mv ${file}_*R2*.fastq ${file}_R2.fastq; done 

Then transferred all “R1.fastq” & “R2.fastq” files to the folder to be submitted to SymPortal

2 SymPortal ITS2 analysis

Information on Symportal output documents here

Formatting notes:

  • Cleaned up file called 188_20211214_03_DBV_20211215T030311.profiles.absolute.abund_and_meta.txt to format for phyloseq

2.1 Setup

## Warning: package 'phyloseq' was built under R version 4.0.3
## Loading required package: lattice
## Loading required package: plyr
## Warning: package 'plyr' was built under R version 4.0.5
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
## 
##     get_legend
## The following object is masked from 'package:plyr':
## 
##     mutate
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.0.5
## This is vegan 2.6-2
## Warning: package 'colorspace' was built under R version 4.0.5
## 
## Attaching package: 'funfuns'
## The following objects are masked from 'package:colorspace':
## 
##     darken, lighten

2.2 Checking for no read samples

counts <- read.csv('flits_symportal_profile_counts.csv',header=TRUE,row.names=1,check.names=FALSE)
#325 samples, 75 type profiles
#+1 sample that's fake 'UK1-O_000' so I can plot UK1-O below

sort(rowSums(counts))
## LK1-I_1199_srad_2018 UK2-O_1105_srad_2018  UK1-O_000_srad_2017 
##                    0                    0                    1 
## UK1-I_1148_ssid_2018 UK1-I_1141_ssid_2018 LK1-I_1192_srad_2018 
##                  494                  550                  714 
## UK1-O_1173_ssid_2018 LK1-O_1229_srad_2018 UK2-I_1074_srad_2018 
##                  720                  873                 1139 
## LK1-I_1196_srad_2018 UK2-I_1085_ssid_2018 LK1-I_1193_srad_2018 
##                 1148                 1334                 1423 
## UK1-I_1144_ssid_2018 UK2-I_1078_srad_2018 UK1-I_1143_ssid_2018 
##                 1604                 1624                 1706 
## UK2-I_1089_ssid_2018 LK1-I_1197_srad_2018 LK1-I_1191_srad_2018 
##                 1757                 1804                 1806 
## UK3-I_1037_ssid_2016 UK1-I_1140_ssid_2018 UK2-O_1118_ssid_2018 
##                 1913                 1914                 1927 
## UK2-I_1084_ssid_2018 UK2-O_1103_srad_2018 LK1-O_1223_srad_2018 
##                 1958                 2005                 2039 
## LK1-O_1226_srad_2018 LK1-O_1232_ssid_2018 LK1-O_1239_ssid_2018 
##                 2046                 2061                 2114 
## UK1-O_1176_ssid_2018 UK2-I_1076_srad_2018 UK2-O_1115_ssid_2018 
##                 2130                 2235                 2241 
## LK1-O_1236_ssid_2018 LK1-I_1194_srad_2018 LK1-O_1221_srad_2018 
##                 2272                 2273                 2290 
## LK1-O_1233_ssid_2018 UK1-I_1149_ssid_2018 UK3-I_1052_srad_2016 
##                 2299                 2324                 2534 
##  LK1-O_928_srad_2017 UK1-I_1142_ssid_2018 UK3-I_1050_srad_2016 
##                 2567                 2570                 2621 
## UK2-O_1116_ssid_2018 LK1-I_1190_srad_2018 UK1-O_1174_ssid_2018 
##                 2685                 2755                 2764 
## UK1-O_1172_ssid_2018 LK1-I_1200_ssid_2018 LK1-O_1235_ssid_2018 
##                 2822                 2830                 2893 
## LK1-I_1195_srad_2018 LK1-I_1207_ssid_2018    LK1-O_8_ssid_2015 
##                 2901                 2981                 3129 
## LK1-I_1209_ssid_2018 UK1-I_1147_ssid_2018   UK2-O_97_ssid_2015 
##                 3134                 3223                 3306 
## UK2-I_1075_srad_2018 LK1-O_1227_srad_2018 UK2-I_1086_ssid_2018 
##                 3350                 3379                 3382 
## UK2-I_1088_ssid_2018  LK1-I_988_srad_2016 UK2-O_1108_srad_2018 
##                 3492                 3611                 3684 
## UK1-I_1145_ssid_2018 UK1-O_1171_ssid_2018 UK2-O_1101_srad_2018 
##                 3815                 3876                 3907 
## UK2-O_1117_ssid_2018 LK1-O_1237_ssid_2018 UK2-O_1119_ssid_2018 
##                 3911                 3918                 3995 
## UK1-O_1177_ssid_2018   LK1-O_27_srad_2015 LK1-I_1204_ssid_2018 
##                 4012                 4066                 4082 
##  UK2-I_846_srad_2017 UK2-I_1072_srad_2018  LK1-I_990_srad_2016 
##                 4114                 4334                 4359 
## UK1-O_1170_ssid_2018 UK3-I_1051_srad_2016 LK1-I_1203_ssid_2018 
##                 4499                 4525                 4627 
## UK2-I_1083_ssid_2018 UK2-O_1107_srad_2018  LK1-O_931_srad_2017 
##                 4675                 4746                 4793 
## LK1-I_1208_ssid_2018    LK1-O_6_ssid_2015 UK2-I_1071_srad_2018 
##                 4795                 4972                 4983 
##   UK1-I_84_srad_2015 LK1-O_1231_ssid_2018 LK1-I_1202_ssid_2018 
##                 5039                 5058                 5148 
## LK1-O_1230_ssid_2018  LK1-I_982_ssid_2016 UK2-I_1082_ssid_2018 
##                 5190                 5315                 5386 
## UK2-I_1081_ssid_2018 LK1-I_1198_srad_2018  UK1-O_783_ssid_2017 
##                 5409                 5417                 5432 
## UK2-O_1113_ssid_2018  UK2-O_800_ssid_2017 LK1-I_1201_ssid_2018 
##                 5434                 5466                 5471 
##  UK1-I_705_ssid_2017 LK1-I_1206_ssid_2018 UK3-I_1054_srad_2016 
##                 5487                 5731                 5836 
## UK3-I_1049_srad_2016 LK1-O_1238_ssid_2018  UK2-O_889_ssid_2017 
##                 5848                 5897                 5907 
## UK1-O_1178_ssid_2018 UK2-O_1109_srad_2018   LK1-O_10_ssid_2015 
##                 6026                 6356                 6411 
##   UK2-O_96_ssid_2015  UK1-O_772_ssid_2017 LK1-O_1222_srad_2018 
##                 6430                 6480                 6560 
##  LK1-O_922_ssid_2017  UK1-O_765_ssid_2017 UK2-O_1112_ssid_2018 
##                 6589                 6621                 6844 
##   LK1-O_26_srad_2015  LK1-O_925_ssid_2017 UK3-O_1026_srad_2016 
##                 7140                 7169                 7210 
## LK1-O_1228_srad_2018   LK1-O_23_srad_2015  UK1-O_768_ssid_2017 
##                 7326                 7359                 7573 
##  UK2-O_882_ssid_2017 UK2-I_1077_srad_2018  LK1-I_965_ssid_2017 
##                 7581                 7655                 7738 
## UK2-I_1087_ssid_2018  UK2-I_853_srad_2017  UK2-O_883_ssid_2017 
##                 7916                 7999                 8120 
##  UK2-I_143_srad_2015   LK1-O_24_srad_2015  UK2-I_147_srad_2015 
##                 8205                 8243                 8319 
##  UK2-I_849_srad_2017  UK2-I_852_srad_2017  UK2-O_113_srad_2015 
##                 8445                 8580                 8619 
##  UK2-O_116_srad_2015   LK1-O_30_srad_2015    LK1-O_2_ssid_2015 
##                 8658                 8674                 9011 
##   UK1-O_35_ssid_2015 UK3-I_1039_ssid_2016  LK1-I_981_ssid_2016 
##                 9045                 9081                 9090 
## UK1-O_1175_ssid_2018  UK1-O_771_ssid_2017  UK2-O_907_srad_2017 
##                 9100                 9127                 9153 
##    LK1-O_4_ssid_2015 UK2-O_1100_srad_2018 UK3-I_1045_ssid_2016 
##                 9193                 9213                 9300 
##  UK2-I_854_srad_2017    LK1-O_3_ssid_2015  UK1-O_763_ssid_2017 
##                 9571                 9707                 9798 
##  UK1-O_769_ssid_2017  UK2-I_850_srad_2017  LK1-I_985_ssid_2016 
##                 9817                 9852                 9966 
##  UK1-I_706_ssid_2017  UK1-I_726_srad_2017 UK2-I_1080_ssid_2018 
##                10038                10091                10142 
##  LK1-I_964_srad_2017  UK1-I_725_ssid_2017  UK2-O_881_ssid_2017 
##                10466                10544                10547 
## UK3-O_1020_srad_2016  UK2-I_130_ssid_2015 UK2-O_1102_srad_2018 
##                10571                10603                10623 
##  UK2-I_148_srad_2015  UK2-I_825_ssid_2017  UK1-O_762_ssid_2017 
##                10632                10812                10875 
## UK3-I_1053_srad_2016  UK1-I_707_ssid_2017  UK1-I_731_srad_2017 
##                10943                10962                10969 
## UK3-I_1048_srad_2016 UK1-O_1179_ssid_2018  LK1-I_978_ssid_2016 
##                11057                11267                11647 
##   UK2-O_92_ssid_2015 UK3-O_1010_ssid_2016 UK1-I_1146_ssid_2018 
##                11696                11701                11838 
##  LK1-I_994_srad_2016  UK2-I_823_ssid_2017   UK1-I_66_ssid_2015 
##                11892                12025                12101 
##    LK1-O_1_ssid_2015   UK1-O_51_srad_2015  LK1-O_929_srad_2017 
##                12224                12328                12360 
##  UK1-I_708_ssid_2017  UK2-O_887_ssid_2017  LK1-I_980_ssid_2016 
##                12601                12832                13180 
##    LK1-O_7_ssid_2015  UK2-I_122_ssid_2015   UK1-I_85_srad_2015 
##                13185                13270                13293 
## UK3-I_1042_ssid_2016  UK1-I_727_srad_2017  UK2-O_119_srad_2015 
##                13311                13317                13392 
##  UK1-I_701_ssid_2017  UK1-I_733_srad_2017 UK3-I_1040_ssid_2016 
##                13398                13429                13466 
##  LK1-I_986_ssid_2016  UK2-I_851_srad_2017 UK2-O_1106_srad_2018 
##                13467                13495                13517 
##  UK2-I_828_ssid_2017   UK1-O_56_srad_2015  UK2-I_808_ssid_2017 
##                13541                13657                13923 
##   UK1-O_40_ssid_2015  UK2-O_886_ssid_2017  UK2-O_112_srad_2015 
##                14077                14106                14293 
##  LK1-I_949_ssid_2017   UK1-O_33_ssid_2015 UK3-O_1008_ssid_2016 
##                14315                14320                14409 
##  LK1-I_962_srad_2017  LK1-O_943_ssid_2017   UK1-O_38_ssid_2015 
##                14467                14502                14652 
## UK3-O_1025_srad_2016  UK1-I_732_srad_2017   UK1-I_81_srad_2015 
##                14871                15047                15071 
##  LK1-O_935_srad_2017  UK2-O_885_ssid_2017  UK1-I_704_ssid_2017 
##                15076                15391                15520 
##  UK2-O_909_srad_2017  UK2-I_149_srad_2015 UK3-O_1022_srad_2016 
##                15831                15856                15887 
##  LK1-O_921_ssid_2017  LK1-O_932_srad_2017  LK1-I_977_ssid_2016 
##                16000                16277                16476 
##  LK1-I_993_srad_2016  LK1-I_979_ssid_2016  LK1-O_917_ssid_2017 
##                16639                16694                16801 
##  UK2-O_888_ssid_2017  UK2-O_111_srad_2015  LK1-I_961_srad_2017 
##                16802                16832                16915 
##  UK2-O_120_srad_2015  UK1-O_766_ssid_2017   UK1-I_69_ssid_2015 
##                17019                17123                17803 
##   UK2-O_98_ssid_2015   UK1-O_31_ssid_2015   UK1-O_32_ssid_2015 
##                17812                17822                17907 
##  UK2-I_121_ssid_2015   UK1-I_86_srad_2015  UK1-O_764_ssid_2017 
##                17913                18157                18208 
##  UK1-O_770_ssid_2017  UK1-O_761_ssid_2017  UK2-I_126_ssid_2015 
##                18255                18471                18513 
##  LK1-O_923_ssid_2017  LK1-O_926_srad_2017 UK3-O_1019_srad_2016 
##                18540                18626                19010 
##   LK1-O_22_srad_2015   UK2-O_95_ssid_2015  UK2-I_128_ssid_2015 
##                19078                19226                19236 
## UK3-O_1018_srad_2016  UK2-I_125_ssid_2015  LK1-I_950_ssid_2017 
##                19258                19595                19684 
##   UK1-O_58_srad_2015   UK2-O_93_ssid_2015  UK1-I_703_ssid_2017 
##                19760                20054                20152 
##  UK2-I_826_ssid_2017  UK2-I_145_srad_2015  UK1-I_709_ssid_2017 
##                20261                20332                20403 
##  LK1-I_948_ssid_2017   UK1-I_65_ssid_2015  LK1-I_984_ssid_2016 
##                20408                20800                20846 
##   UK1-O_36_ssid_2015  LK1-O_919_ssid_2017  UK2-I_824_ssid_2017 
##                21174                21265                21277 
##   UK1-I_68_ssid_2015  LK1-I_951_ssid_2017  UK2-I_123_ssid_2015 
##                21408                21520                21552 
##  UK2-O_915_srad_2017  UK1-I_711_ssid_2017   UK2-O_91_ssid_2015 
##                21573                21663                21798 
##   UK2-O_94_ssid_2015  UK2-I_127_ssid_2015    LK1-O_5_ssid_2015 
##                21872                22036                22255 
##  UK2-I_809_ssid_2017  UK2-O_100_ssid_2015 UK3-I_1041_ssid_2016 
##                22499                22568                22603 
##  LK1-I_963_srad_2017  UK2-I_829_ssid_2017  LK1-O_916_ssid_2017 
##                22711                22902                22995 
##  LK1-O_927_srad_2017 UK3-O_1023_srad_2016   UK1-I_62_ssid_2015 
##                23150                23173                23633 
##   UK1-O_54_srad_2015  UK2-I_129_ssid_2015  LK1-I_976_srad_2017 
##                24264                24434                25940 
##  UK1-I_735_srad_2017  UK2-O_884_ssid_2017  LK1-I_956_srad_2017 
##                26183                26322                26615 
##   UK1-O_37_ssid_2015  UK2-I_124_ssid_2015   UK1-I_82_srad_2015 
##                26651                27593                27771 
##  UK2-O_118_srad_2015  UK1-O_767_ssid_2017  LK1-I_958_srad_2017 
##                28102                28250                28449 
##  UK2-I_855_srad_2017 UK3-I_1046_ssid_2016 UK3-O_1017_srad_2016 
##                28604                28623                29057 
##  UK2-O_913_srad_2017   UK1-O_57_srad_2015 UK3-O_1014_ssid_2016 
##                29116                29330                29454 
## UK3-O_1011_ssid_2016  LK1-I_947_ssid_2017  UK2-O_890_ssid_2017 
##                30004                30367                30503 
##  UK2-I_807_ssid_2017   UK1-I_64_ssid_2015   UK1-O_39_ssid_2015 
##                31352                31518                32161 
## UK3-O_1013_ssid_2016  LK1-I_955_ssid_2017  UK2-O_114_srad_2015 
##                32219                32224                32661 
##  UK2-O_798_ssid_2017  UK2-I_827_ssid_2017   UK1-I_63_ssid_2015 
##                32936                33128                34025 
##  UK1-I_702_ssid_2017 UK3-O_1016_ssid_2016  UK2-O_115_srad_2015 
##                34417                35721                35763 
##   LK1-O_21_srad_2015  UK2-O_914_srad_2017   UK1-O_34_ssid_2015 
##                36207                37720                38141 
##  UK2-O_912_srad_2017 UK3-O_1007_ssid_2016  UK2-O_117_srad_2015 
##                38907                39109                39575 
##   UK1-I_70_ssid_2015   UK1-I_61_ssid_2015   UK1-I_67_ssid_2015 
##                39625                39692                39935 
##  UK2-I_821_ssid_2017 UK3-O_1024_srad_2016  UK1-I_710_ssid_2017 
##                40254                41155                41230 
##  UK2-I_822_ssid_2017  LK1-I_983_ssid_2016  UK2-O_908_srad_2017 
##                41969                43004                43349 
##  LK1-I_996_srad_2016 UK3-O_1009_ssid_2016  UK2-I_847_srad_2017 
##                44107                44242                44413 
##  LK1-I_992_srad_2016   UK1-I_89_srad_2015  UK2-O_820_ssid_2017 
##                44988                45207                45446 
##  UK2-I_150_srad_2015 UK3-O_1015_ssid_2016  LK1-I_995_srad_2016 
##                45889                46857                47563 
##  UK2-I_848_srad_2017   UK2-O_99_ssid_2015  UK2-I_830_ssid_2017 
##                49830                50651                69634 
##  UK2-I_144_srad_2015  UK2-I_801_ssid_2017 
##                77097                88938
counts.no0 <- counts[rowSums(counts)>0,]
#2 that are 0, removing them now: LK1-I_1199_srad_2018, UK2-O_1105_srad_2018

2.3 Phyloseq object - type profiles

Ran once then saved to re-read in later

# import dataframe holding sample information
samdf <- read.csv("fl_sample_info - ITS2_all.csv")
head(samdf)
##   sample_num full_site site zone site_zone  sample  spp year sequencing time
## 1          1   E Sambo  LK1    O     LK1-O LK1-O_1 ssid 2015       ITS2  Pre
## 2          2   E Sambo  LK1    O     LK1-O LK1-O_2 ssid 2015       ITS2  Pre
## 3          3   E Sambo  LK1    O     LK1-O LK1-O_3 ssid 2015       ITS2  Pre
## 4          4   E Sambo  LK1    O     LK1-O LK1-O_4 ssid 2015       ITS2  Pre
## 5          5   E Sambo  LK1    O     LK1-O LK1-O_5 ssid 2015       ITS2  Pre
## 6          6   E Sambo  LK1    O     LK1-O LK1-O_6 ssid 2015       ITS2  Pre
##   season
## 1 Spring
## 2 Spring
## 3 Spring
## 4 Spring
## 5 Spring
## 6 Spring
samdf$full_sample <- paste0(samdf$sample,"_",samdf$spp,"_",samdf$year)
rownames(samdf) <- samdf$full_sample

# import taxa info
taxa <- read.csv("flits_symportal_taxa.csv",header=TRUE,row.names=1)
#rownames(taxa) <- as.factor(taxa$ITS2_type_profile)
mtaxa <- as.matrix(taxa)

# import counts (absolute abundance from its2 type profiles)
mcounts <- as.matrix(counts.no0)

# Construct phyloseq object 
ps <- phyloseq(otu_table(mcounts, taxa_are_rows=FALSE),
               sample_data(samdf),
               tax_table(mtaxa))

ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 75 taxa and 324 samples ]
## sample_data() Sample Data:       [ 324 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 75 taxa by 5 taxonomic ranks ]
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 75 taxa and 324 samples ]
# sample_data() Sample Data:       [ 324 samples by 12 sample variables ]
# tax_table()   Taxonomy Table:    [ 75 taxa by 5 taxonomic ranks ]

#saveRDS(ps,file="ps.its2.rds")

2.4 Re-read in data

ps <- readRDS("ps.its2.RDS")

2.5 Bar plots by majority sequence

2.5.1 Assessing type profile distances

Note - used file “20211215T030311_unifrac_profiles_PCoA_coords_B_sqrt.csv”, made the name shorter. I also removed the ‘proportion explained’ line at the end. Also added the first type that shows up as its own column, ‘first_type’

Proportion explained: PC1 0.584828597, PC2 0.146734687 For bray curtis: PC1 0.231, PC2 0.139

library(ggrepel)

setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")

uni.cord <- read.csv("unifrac_profiles_PCoA_coords_B_sqrt.csv")

ggplot(uni.cord,aes(x=PC1,y=PC2,color=type_collapsed,label=sample))+
  geom_point()+
  xlab("PC1 (58.5%)")+
  ylab("PC2 (14.7%)")+
  #geom_text_repel(max.overlaps=20)+
  stat_ellipse()+
  theme_bw()
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Warning: Removed 3 row(s) containing missing values (geom_path).

#separated B5-1* & B5-2* by the 0 on PC2, except for 1 that was just below 0

#ggsave(file="pca_grouped.pdf",w=12,h=10)

#checking bray curtis
bray.cord <- read.csv("braycurtis_profiles_PCoA_coords_B_sqrt.csv")

ggplot(bray.cord,aes(x=PC1,y=PC2,color=type_bray,label=sample))+
  geom_point()+
  #geom_text_repel(max.overlaps=20)+
  stat_ellipse()
## Too few points to calculate an ellipse
## Warning: Removed 1 row(s) containing missing values (geom_path).

Important setup

ps.maj <- tax_glom(ps,"Majority_ITS2_collapsed")
#ps.maj <- tax_glom(ps,"Majority_ITS2_sequence_bray")

ntaxa(ps); ntaxa(ps.maj) #18 for unifrac sids, 14 for rads
## [1] 75
## [1] 14
#16 for bray sids

ps.maj <- subset_samples(ps.maj,site!="UK3")

#siderea
ps.maj.sid <- subset_samples(ps.maj,spp=="ssid")

ps.maj.sid.pre <- subset_samples(ps.maj.sid,time=="Pre")
ps.maj.sid.pre.lk1 <- subset_samples(ps.maj.sid.pre,site=="LK1")
ps.maj.sid.pre.uk1 <- subset_samples(ps.maj.sid.pre,site=="UK1")
ps.maj.sid.pre.uk2 <- subset_samples(ps.maj.sid.pre,site=="UK2")

ps.maj.sid.irm <- subset_samples(ps.maj.sid,time=="Irma")
ps.maj.sid.irm.lk1 <- subset_samples(ps.maj.sid.irm,site=="LK1")
ps.maj.sid.irm.uk1 <- subset_samples(ps.maj.sid.irm,site=="UK1")
ps.maj.sid.irm.uk2 <- subset_samples(ps.maj.sid.irm,site=="UK2")

ps.maj.sid.pos <- subset_samples(ps.maj.sid,time=="Post")
ps.maj.sid.pos.lk1 <- subset_samples(ps.maj.sid.pos,site=="LK1")
ps.maj.sid.pos.uk1 <- subset_samples(ps.maj.sid.pos,site=="UK1")
ps.maj.sid.pos.uk2 <- subset_samples(ps.maj.sid.pos,site=="UK2")

#radians
ps.maj.rad <- subset_samples(ps.maj,spp=="srad")

ps.maj.rad.pre <- subset_samples(ps.maj.rad,time=="Pre")
ps.maj.rad.pre.lk1 <- subset_samples(ps.maj.rad.pre,site=="LK1")
ps.maj.rad.pre.uk1 <- subset_samples(ps.maj.rad.pre,site=="UK1")
ps.maj.rad.pre.uk2 <- subset_samples(ps.maj.rad.pre,site=="UK2")

ps.maj.rad.irm <- subset_samples(ps.maj.rad,time=="Irma")
ps.maj.rad.irm.lk1 <- subset_samples(ps.maj.rad.irm,site=="LK1")
ps.maj.rad.irm.uk1 <- subset_samples(ps.maj.rad.irm,site=="UK1")
ps.maj.rad.irm.uk2 <- subset_samples(ps.maj.rad.irm,site=="UK2")

ps.maj.rad.pos <- subset_samples(ps.maj.rad,time=="Post")
ps.maj.rad.pos.lk1 <- subset_samples(ps.maj.rad.pos,site=="LK1")
#ps.maj.rad.pos.uk1 <- subset_samples(ps.maj.rad.pos,site=="UK1")
#not available
ps.maj.rad.pos.uk2 <- subset_samples(ps.maj.rad.pos,site=="UK2")

Removing black bars from phyloseq plot

plot_bar2 <-  function (physeq, x = "Sample", y = "Abundance", fill = NULL, title = NULL, facet_grid = NULL, border_color = NA) 
{
  mdf = psmelt(physeq)
  p = ggplot(mdf, aes_string(x = x, y = y, fill = fill))
  p = p + geom_bar(stat = "identity", position = "stack",  color = border_color)
  p = p + theme(axis.text.x = element_text(angle = -90, hjust = 0))
  if (!is.null(facet_grid)) {
    p <- p + facet_grid(facet_grid)
  }
  if (!is.null(title)) {
    p <- p + ggtitle(title)
  }
  return(p)
}

maj.seq.colors <- c("#474747","#7A7A7A","#A8A8A8","#00294D","#195481","#5882AF","#9AB0CE","#D6DDE9","#D3E0D2","#95B792","#548C4E","#1F5D12","#002F00","#B8773A")

2.5.2 The plots

2.5.2.1 S. siderea

LOWER KEYS 1

ps.maj.sid.pre.lk1.rel <- transform_sample_counts(ps.maj.sid.pre.lk1, function(x) x / sum(x))
ps.maj.sid.pre.lk1.rel.zone <- merge_samples(ps.maj.sid.pre.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pre.lk1.zone <- transform_sample_counts(ps.maj.sid.pre.lk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.pre.lk1.zone <- plot_bar2(ps.maj.sid.pre.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pre.lk1.zone

ps.maj.sid.irm.lk1.rel <- transform_sample_counts(ps.maj.sid.irm.lk1, function(x) x / sum(x))
ps.maj.sid.irm.lk1.rel.zone <- merge_samples(ps.maj.sid.irm.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.irm.lk1.zone <- transform_sample_counts(ps.maj.sid.irm.lk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.irm.lk1.zone <- plot_bar2(ps.maj.sid.irm.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.irm.lk1.zone

ps.maj.sid.pos.lk1.rel <- transform_sample_counts(ps.maj.sid.pos.lk1, function(x) x / sum(x))
ps.maj.sid.pos.lk1.rel.zone <- merge_samples(ps.maj.sid.pos.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pos.lk1.zone <- transform_sample_counts(ps.maj.sid.pos.lk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.pos.lk1.zone <- plot_bar2(ps.maj.sid.pos.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pos.lk1.zone

#ggarrange(gg.maj.sid.pre.lk1.zone,gg.maj.sid.irm.lk1.zone,gg.maj.sid.pos.lk1.zone,nrow=1)

UPPER KEYS 1

ps.maj.sid.pre.uk1.rel <- transform_sample_counts(ps.maj.sid.pre.uk1, function(x) x / sum(x))
ps.maj.sid.pre.uk1.rel.zone <- merge_samples(ps.maj.sid.pre.uk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pre.uk1.zone <- transform_sample_counts(ps.maj.sid.pre.uk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.pre.uk1.zone <- plot_bar2(ps.maj.sid.pre.uk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pre.uk1.zone

ps.maj.sid.irm.uk1.rel <- transform_sample_counts(ps.maj.sid.irm.uk1, function(x) x / sum(x))
ps.maj.sid.irm.uk1.rel.zone <- merge_samples(ps.maj.sid.irm.uk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.irm.uk1.zone <- transform_sample_counts(ps.maj.sid.irm.uk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.irm.uk1.zone <- plot_bar2(ps.maj.sid.irm.uk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.irm.uk1.zone

ps.maj.sid.pos.uk1.rel <- transform_sample_counts(ps.maj.sid.pos.uk1, function(x) x / sum(x))
ps.maj.sid.pos.uk1.rel.zone <- merge_samples(ps.maj.sid.pos.uk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pos.uk1.zone <- transform_sample_counts(ps.maj.sid.pos.uk1.rel.zone, function(x) x / sum(x))

gg.maj.sid.pos.uk1.zone <- plot_bar2(ps.maj.sid.pos.uk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pos.uk1.zone

#ggarrange(gg.maj.sid.pre.uk1.zone,gg.maj.sid.irm.uk1.zone,gg.maj.sid.pos.uk1.zone,nrow=1)

UPPER KEYS 2

ps.maj.sid.pre.uk2.rel <- transform_sample_counts(ps.maj.sid.pre.uk2, function(x) x / sum(x))
ps.maj.sid.pre.uk2.rel.zone <- merge_samples(ps.maj.sid.pre.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pre.uk2.zone <- transform_sample_counts(ps.maj.sid.pre.uk2.rel.zone, function(x) x / sum(x))

gg.maj.sid.pre.uk2.zone <- plot_bar2(ps.maj.sid.pre.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pre.uk2.zone

ps.maj.sid.irm.uk2.rel <- transform_sample_counts(ps.maj.sid.irm.uk2, function(x) x / sum(x))
ps.maj.sid.irm.uk2.rel.zone <- merge_samples(ps.maj.sid.irm.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.irm.uk2.zone <- transform_sample_counts(ps.maj.sid.irm.uk2.rel.zone, function(x) x / sum(x))

gg.maj.sid.irm.uk2.zone <- plot_bar2(ps.maj.sid.irm.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.irm.uk2.zone

ps.maj.sid.pos.uk2.rel <- transform_sample_counts(ps.maj.sid.pos.uk2, function(x) x / sum(x))
ps.maj.sid.pos.uk2.rel.zone <- merge_samples(ps.maj.sid.pos.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pos.uk2.zone <- transform_sample_counts(ps.maj.sid.pos.uk2.rel.zone, function(x) x / sum(x))

gg.maj.sid.pos.uk2.zone <- plot_bar2(ps.maj.sid.pos.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.sid.pos.uk2.zone

#ggarrange(gg.maj.sid.pre.uk2.zone,gg.maj.sid.irm.uk2.zone,gg.maj.sid.pos.uk2.zone,nrow=1)

All the sites

ggarrange(gg.maj.sid.pre.lk1.zone,gg.maj.sid.irm.lk1.zone,gg.maj.sid.pos.lk1.zone,gg.maj.sid.pre.uk1.zone,gg.maj.sid.irm.uk1.zone,gg.maj.sid.pos.uk1.zone,gg.maj.sid.pre.uk2.zone,gg.maj.sid.irm.uk2.zone,gg.maj.sid.pos.uk2.zone,nrow=3,ncol=3,common.legend=T,legend="right")

#ggsave("gg.its2.sids.pdf",h=8,w=7)

2.5.2.2 S. radians

LOWER KEYS 1

ps.maj.rad.pre.lk1.rel <- transform_sample_counts(ps.maj.rad.pre.lk1, function(x) x / sum(x))
ps.maj.rad.pre.lk1.rel.zone <- merge_samples(ps.maj.rad.pre.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.pre.lk1.zone <- transform_sample_counts(ps.maj.rad.pre.lk1.rel.zone, function(x) x / sum(x))

gg.maj.rad.pre.lk1.zone <- plot_bar2(ps.maj.rad.pre.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.pre.lk1.zone

ps.maj.rad.irm.lk1.rel <- transform_sample_counts(ps.maj.rad.irm.lk1, function(x) x / sum(x))
ps.maj.rad.irm.lk1.rel.zone <- merge_samples(ps.maj.rad.irm.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.irm.lk1.zone <- transform_sample_counts(ps.maj.rad.irm.lk1.rel.zone, function(x) x / sum(x))

gg.maj.rad.irm.lk1.zone <- plot_bar2(ps.maj.rad.irm.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.irm.lk1.zone

ps.maj.rad.pos.lk1.rel <- transform_sample_counts(ps.maj.rad.pos.lk1, function(x) x / sum(x))
ps.maj.rad.pos.lk1.rel.zone <- merge_samples(ps.maj.rad.pos.lk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.pos.lk1.zone <- transform_sample_counts(ps.maj.rad.pos.lk1.rel.zone, function(x) x / sum(x))

gg.maj.rad.pos.lk1.zone <- plot_bar2(ps.maj.rad.pos.lk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.pos.lk1.zone

#ggarrange(gg.maj.rad.pre.lk1.zone,gg.maj.rad.irm.lk1.zone,gg.maj.rad.pos.lk1.zone,nrow=1)

UPPER KEYS 1

ps.maj.rad.pre.uk1.rel <- transform_sample_counts(ps.maj.rad.pre.uk1, function(x) x / sum(x))
ps.maj.rad.pre.uk1.rel.zone <- merge_samples(ps.maj.rad.pre.uk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.pre.uk1.zone <- transform_sample_counts(ps.maj.rad.pre.uk1.rel.zone, function(x) x / sum(x))

gg.maj.rad.pre.uk1.zone <- plot_bar2(ps.maj.rad.pre.uk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.pre.uk1.zone

ps.maj.rad.irm.uk1.rel <- transform_sample_counts(ps.maj.rad.irm.uk1, function(x) x / sum(x))
ps.maj.rad.irm.uk1.rel.zone <- merge_samples(ps.maj.rad.irm.uk1.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.irm.uk1.zone <- transform_sample_counts(ps.maj.rad.irm.uk1.rel.zone, function(x) x / sum(x))

gg.maj.rad.irm.uk1.zone <- plot_bar2(ps.maj.rad.irm.uk1.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.irm.uk1.zone

#uk1 not available for Post

#ggarrange(gg.maj.rad.pre.uk1.zone,gg.maj.rad.irm.uk1.zone,gg.maj.rad.pos.uk1.zone,nrow=1)

UPPER KEYS 2

ps.maj.rad.pre.uk2.rel <- transform_sample_counts(ps.maj.rad.pre.uk2, function(x) x / sum(x))
ps.maj.rad.pre.uk2.rel.zone <- merge_samples(ps.maj.rad.pre.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.pre.uk2.zone <- transform_sample_counts(ps.maj.rad.pre.uk2.rel.zone, function(x) x / sum(x))

gg.maj.rad.pre.uk2.zone <- plot_bar2(ps.maj.rad.pre.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.pre.uk2.zone

ps.maj.rad.irm.uk2.rel <- transform_sample_counts(ps.maj.rad.irm.uk2, function(x) x / sum(x))
ps.maj.rad.irm.uk2.rel.zone <- merge_samples(ps.maj.rad.irm.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.irm.uk2.zone <- transform_sample_counts(ps.maj.rad.irm.uk2.rel.zone, function(x) x / sum(x))

gg.maj.rad.irm.uk2.zone <- plot_bar2(ps.maj.rad.irm.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturb.")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.irm.uk2.zone

ps.maj.rad.pos.uk2.rel <- transform_sample_counts(ps.maj.rad.pos.uk2, function(x) x / sum(x))
ps.maj.rad.pos.uk2.rel.zone <- merge_samples(ps.maj.rad.pos.uk2.rel, "zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.rad.pos.uk2.zone <- transform_sample_counts(ps.maj.rad.pos.uk2.rel.zone, function(x) x / sum(x))

gg.maj.rad.pos.uk2.zone <- plot_bar2(ps.maj.rad.pos.uk2.zone, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_x_discrete(labels=c("Inshore","Offshore"))+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")+
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  theme(axis.text.x=element_text(angle=45))
  
gg.maj.rad.pos.uk2.zone

#ggarrange(gg.maj.rad.pre.uk2.zone,gg.maj.rad.irm.uk2.zone,gg.maj.rad.pos.uk2.zone,nrow=1)

All the sites

ggarrange(gg.maj.rad.pre.lk1.zone,gg.maj.rad.irm.lk1.zone,gg.maj.rad.pos.lk1.zone,gg.maj.rad.pre.uk1.zone,gg.maj.rad.irm.uk1.zone,gg.maj.rad.pre.uk2.zone,gg.maj.rad.irm.uk2.zone,gg.maj.rad.pos.uk2.zone,nrow=3,ncol=3,common.legend=T,legend="right")

#ggsave("gg.its2.rads.pdf",h=8,w=7)

2.5.3 Different version of the same type profile plot

ps.maj.sid.pre.rel <- transform_sample_counts(ps.maj.sid.pre, function(x) x / sum(x))
ps.maj.sid.pre.rel.sz <- merge_samples(ps.maj.sid.pre.rel, "site_zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pre.rel.sz <- transform_sample_counts(ps.maj.sid.pre.rel.sz, function(x) x / sum(x))

gg.bar.pre <- plot_bar2(ps.maj.sid.pre.rel.sz, fill="Majority_ITS2_collapsed")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_fill_manual(values=maj.seq.colors)+
  geom_bar(stat = "identity")
gg.bar.pre

ps.maj.sid.irm.rel <- transform_sample_counts(ps.maj.sid.irm, function(x) x / sum(x))
ps.maj.sid.irm.rel.sz <- merge_samples(ps.maj.sid.irm.rel, "site_zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.irm.rel.sz <- transform_sample_counts(ps.maj.sid.irm.rel.sz, function(x) x / sum(x))

gg.bar.irm <- plot_bar(ps.maj.sid.irm.rel.sz, fill="Majority_ITS2_collapsed")+
    ggtitle("Disturbance")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
    scale_fill_manual(values=maj.seq.colors)
gg.bar.irm

ps.maj.sid.pos.rel <- transform_sample_counts(ps.maj.sid.pos, function(x) x / sum(x))
ps.maj.sid.pos.rel.sz <- merge_samples(ps.maj.sid.pos.rel, "site_zone")
## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion

## Warning in asMethod(object): NAs introduced by coercion
ps.maj.sid.pos.rel.sz <- transform_sample_counts(ps.maj.sid.pos.rel.sz, function(x) x / sum(x))

gg.bar.pos <- plot_bar(ps.maj.sid.pos.rel.sz, fill="Majority_ITS2_collapsed")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
    scale_fill_manual(values=maj.seq.colors)
gg.bar.pos

gg.bar.its2 <- ggarrange(gg.bar.pre,gg.bar.irm,gg.bar.pos,ncol=3,common.legend=T,legend="right")
gg.bar.its2

#ggsave("bar.its2.sids.pdf",width=9,height=4)

Type profile - Bray-Curtis distance

maj.seq.colors <- c("#1B1B1B","#969696","#F9F9F9","#00294D","#195481","#5882AF","#9AB0CE","#D6DDE9","#D3E0D2","#95B792","#548C4E","#1F5D12","#002F00","#B8773A")
                         
ps.maj.sid.pre.rel <- transform_sample_counts(ps.maj.sid.pre, function(x) x / sum(x))
ps.maj.sid.pre.rel.sz <- merge_samples(ps.maj.sid.pre.rel, "site_zone")
ps.maj.sid.pre.rel.sz <- transform_sample_counts(ps.maj.sid.pre.rel.sz, function(x) x / sum(x))

gg.bar.pre <- plot_bar(ps.maj.sid.pre.rel.sz, fill="Majority_ITS2_sequence_bray")+
    ggtitle("Pre")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
  scale_fill_manual(values=maj.seq.colors)
gg.bar.pre

ps.maj.sid.irm.rel <- transform_sample_counts(ps.maj.sid.irm, function(x) x / sum(x))
ps.maj.sid.irm.rel.sz <- merge_samples(ps.maj.sid.irm.rel, "site_zone")
ps.maj.sid.irm.rel.sz <- transform_sample_counts(ps.maj.sid.irm.rel.sz, function(x) x / sum(x))

gg.bar.irm <- plot_bar(ps.maj.sid.irm.rel.sz, fill="Majority_ITS2_sequence_bray")+
    ggtitle("Disturbance")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
    scale_fill_manual(values=maj.seq.colors)
gg.bar.irm

ps.maj.sid.pos.rel <- transform_sample_counts(ps.maj.sid.pos, function(x) x / sum(x))
ps.maj.sid.pos.rel.sz <- merge_samples(ps.maj.sid.pos.rel, "site_zone")
ps.maj.sid.pos.rel.sz <- transform_sample_counts(ps.maj.sid.pos.rel.sz, function(x) x / sum(x))

gg.bar.pos <- plot_bar(ps.maj.sid.pos.rel.sz, fill="Majority_ITS2_sequence_bray")+
    ggtitle("Post")+
    theme_cowplot()+
    theme(axis.text.x=element_text(angle=90,hjust=1))+
    xlab("")+
    scale_fill_manual(values=maj.seq.colors)
gg.bar.pos

gg.bar.its2 <- ggarrange(gg.bar.pre,gg.bar.irm,gg.bar.pos,ncol=3,common.legend=T,legend="right")
gg.bar.its2
#ggsave("bar.its2.sids.pdf",width=9,height=4)

2.6 Bar plot by type profile

2.6.1 Setup - packages

#remotes::install_github("KarstensLab/microshades")
library("microshades")
#remotes::install_github("mikemc/speedyseq")
library("speedyseq")
## 
## Attaching package: 'speedyseq'
## The following objects are masked from 'package:phyloseq':
## 
##     filter_taxa, plot_bar, plot_heatmap, plot_tree, psmelt, tax_glom,
##     tip_glom, transform_sample_counts
library(ggpubr)
library(cowplot)

#setwd("/Volumes/GoogleDrive-104519233854090018057/My Drive/Flirma/flirma_networks/fl_its2")
setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")

2.6.2 Setup - data

ps <- readRDS("ps.its2.rds")

ps <- subset_samples(ps,site!="UK3")

ps.sid <- subset_samples(ps,spp=="ssid")
ps.rad <- subset_samples(ps,spp=="srad")

ps.sid.no0 <- prune_taxa(taxa_sums(ps.sid)!=0,ps.sid)
ps.rad.no0 <- prune_taxa(taxa_sums(ps.rad)!=0,ps.rad)

2.6.3 Microshades!

2.6.3.1 Microshades plot - sids

# Use microshades function prep_mdf to agglomerate, normalize, and melt the phyloseq object
mdf.ps.sid <- prep_mdf(ps.sid.no0,subgroup_level="Majority_ITS2_collapsed")

type_groups <- c(unique(mdf.ps.sid$Clade))

# Create a color object for the specified data
col.mdf.ps.sid <- create_color_dfs(mdf.ps.sid, top_orientation = FALSE,group_level="Clade",subgroup_level="Majority_ITS2_collapsed",selected_groups=c("A","B","C","D"),cvd=TRUE)

# Extract
col.mdf.ps.sid.m <- col.mdf.ps.sid$mdf
col.mdf.ps.sid.c <- col.mdf.ps.sid$cdf

col.mdf.ps.sid.m$zone <- sub("I","Inshore",col.mdf.ps.sid.m$zone)
col.mdf.ps.sid.m$zone <- sub("O","Offshore",col.mdf.ps.sid.m$zone)

col.mdf.ps.sid.m$time <- factor(col.mdf.ps.sid.m$time,levels=c("Pre","Irma","Post"))

col.mdf.ps.sid.c.new <- color_reassign(col.mdf.ps.sid.c,
                          group_assignment = c("A","B","C","D"),
                          color_assignment = c("micro_gray","micro_cvd_blue","micro_cvd_green","micro_brown"),group_level="Clade")

col.mdf.ps.sid.m.lk1 <- subset(col.mdf.ps.sid.m,site=="LK1")

ms.plot.lk1 <- plot_microshades(col.mdf.ps.sid.m.lk1, col.mdf.ps.sid.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Lower Keys 1")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
  #theme(axis.text.x=element_text(angle=45,hjust=1))
ms.plot.lk1

col.mdf.ps.sid.m.uk1 <- subset(col.mdf.ps.sid.m,site=="UK1")

ms.plot.uk1 <- plot_microshades(col.mdf.ps.sid.m.uk1, col.mdf.ps.sid.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Upper Keys 1")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
ms.plot.uk1

col.mdf.ps.sid.m.uk2 <- subset(col.mdf.ps.sid.m,site=="UK2")

ms.plot.uk2 <- plot_microshades(col.mdf.ps.sid.m.uk2, col.mdf.ps.sid.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Upper Keys 2")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
ms.plot.uk2

ms.plot.sid <- ggarrange(ms.plot.lk1,ms.plot.uk1,ms.plot.uk2,nrow=3,ncol=1,labels="AUTO")
annotate_figure(ms.plot.sid,top=text_grob("S. siderea",face="italic",size=16))

#ggsave("ms.plot.sid.pdf",height=10,width=8)

2.6.3.2 Microshades plot - rads

# Use microshades function prep_mdf to agglomerate, normalize, and melt the phyloseq object
mdf.ps.rad <- prep_mdf(ps.rad.no0,subgroup_level="Majority_ITS2_collapsed")

type_groups <- c(unique(mdf.ps.rad$Clade))

# Create a color object for the specified data
col.mdf.ps.rad <- create_color_dfs(mdf.ps.rad, top_orientation = FALSE,group_level="Clade",subgroup_level="Majority_ITS2_collapsed",selected_groups=c("A","B","C","D"),cvd=TRUE)

# Extract
col.mdf.ps.rad.m <- col.mdf.ps.rad$mdf
col.mdf.ps.rad.c <- col.mdf.ps.rad$cdf

col.mdf.ps.rad.m$zone <- sub("I","Inshore",col.mdf.ps.rad.m$zone)
col.mdf.ps.rad.m$zone <- sub("O","Offshore",col.mdf.ps.rad.m$zone)

col.mdf.ps.rad.m$time <- factor(col.mdf.ps.rad.m$time,levels=c("Pre","Irma","Post"))

col.mdf.ps.rad.c.new <- color_reassign(col.mdf.ps.rad.c,
                          group_assignment = c("A","B","C","D"),
                          color_assignment = c("micro_gray","micro_cvd_blue","micro_cvd_green","micro_brown"),group_level="Clade")

col.mdf.ps.rad.m.lk1 <- subset(col.mdf.ps.rad.m,site=="LK1")

ms.plot.lk1 <- plot_microshades(col.mdf.ps.rad.m.lk1, col.mdf.ps.rad.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Lower Keys 1")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
ms.plot.lk1

col.mdf.ps.rad.m.uk1 <- subset(col.mdf.ps.rad.m,site=="UK1")

ms.plot.uk1 <- plot_microshades(col.mdf.ps.rad.m.uk1, col.mdf.ps.rad.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Upper Keys 1")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
ms.plot.uk1

col.mdf.ps.rad.m.uk2 <- subset(col.mdf.ps.rad.m,site=="UK2")

ms.plot.uk2 <- plot_microshades(col.mdf.ps.rad.m.uk2, col.mdf.ps.rad.c.new)+
  facet_wrap(zone~time,scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x=element_blank())+
  ggtitle("Upper Keys 1")+
  guides(color=guide_legend("Majority ITS2 sequence"),fill=guide_legend("Majority ITS2 sequence")) 
ms.plot.uk2

ms.plot.rad <- ggarrange(ms.plot.lk1,ms.plot.uk1,ms.plot.uk2,nrow=3,ncol=1)
annotate_figure(ms.plot.rad,top=text_grob("S. radians",face="italic",size=16))

#ggsave("ms.plot.rad.pdf",height=10,width=8)

3 PCoAs/stats

3.1 Setup

library(cowplot)
library(phyloseq)
library(vegan)
library(ggplot2)

setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
ps.pm <- readRDS("ps.its2.postmed.rds")

Run once, then skip

counts <- read.csv('flits_postmed.csv',header=TRUE,row.names=1,check.names=FALSE)
#325 samples, 1268 post med seqs
#in relative abundance format

# import dataframe holding sample information
samdf <- read.csv("fl_sample_info - ITS2_all.csv")
head(samdf)
samdf$full_sample <- paste0(samdf$sample,"_",samdf$spp,"_",samdf$year)
rownames(samdf) <- samdf$full_sample

# import counts (relative abundance from its2 type profiles)
mcounts <- as.matrix(counts)

# Construct phyloseq object 
ps.pm1 <- phyloseq(otu_table(mcounts, taxa_are_rows=FALSE),
               sample_data(samdf))

ps.pm1

ps.pm <- subset_samples(ps.pm1,site!="UK3")
ps.pm
# phyloseq-class experiment-level object
# otu_table()   OTU Table:         [ 1268 taxa and 293 samples ]
# sample_data() Sample Data:       [ 293 samples by 12 sample variables ]

#saveRDS(ps.pm,file="ps.its2.postmed.rds")

3.2 Stats by time - sids

3.2.1 Setup

ps.pm.sid <- subset_samples(ps.pm,spp=="ssid")

##for time stats
ps.pm.sid.lk1 <- subset_samples(ps.pm.sid,site=="LK1")
ps.pm.sid.uk1 <- subset_samples(ps.pm.sid,site=="UK1")
ps.pm.sid.uk2 <- subset_samples(ps.pm.sid,site=="UK2")

ps.pm.sid.lk1i <- subset_samples(ps.pm.sid.lk1,zone=="I")
ps.pm.sid.lk1o <- subset_samples(ps.pm.sid.lk1,zone=="O")

ps.pm.sid.uk1i <- subset_samples(ps.pm.sid.uk1,zone=="I")
ps.pm.sid.uk1o <- subset_samples(ps.pm.sid.uk1,zone=="O")

ps.pm.sid.uk2i <- subset_samples(ps.pm.sid.uk2,zone=="I")
ps.pm.sid.uk2o <- subset_samples(ps.pm.sid.uk2,zone=="O")

3.2.2 Time - LK1 sids

##time - inshore
seq.sid.lk1i <- data.frame(ps.pm.sid.lk1i@otu_table)
samdf.sid.lk1i <- data.frame(ps.pm.sid.lk1i@sam_data)

dist.sid.lk1i <- vegdist(seq.sid.lk1i)

bet.sid.lk1i <- betadisper(dist.sid.lk1i,samdf.sid.lk1i$time)
#anova(bet.ps.pm)
permutest(bet.sid.lk1i,pairwise=TRUE,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.08061 0.040304 0.9945    999  0.375
## Residuals 23 0.93215 0.040528                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          Irma     Post   Pre
## Irma          0.029000 0.590
## Post 0.034299          0.446
## Pre  0.615872 0.442412
#post & Irma differ
#          Irma     Post   Pre
# Irma          0.032000 0.603
# Post 0.034299          0.427
# Pre  0.615872 0.442412      
plot(bet.sid.lk1i) 

adonis2(seq.sid.lk1i ~ time, data=samdf.sid.lk1i, permutations=999) #p **
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.lk1i ~ time, data = samdf.sid.lk1i, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)   
## time      2   1.8496 0.25275 3.8898  0.004 **
## Residual 23   5.4683 0.74725                 
## Total    25   7.3179 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.sid.lk1i, factors=samdf.sid.lk1i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs   F.Model         R2 p.value p.adjusted
## 1 Post vs Irma 3.6451374 0.20658028   0.002      0.002
## 2  Post vs Pre 7.7147744 0.31215233   0.001      0.001
## 3  Irma vs Pre 0.7382707 0.04690926   0.441      0.441
##post different than other two
#          pairs   F.Model         R2 p.value p.adjusted
# 1 Post vs Irma 3.6451374 0.20658028   0.004      0.004
# 2  Post vs Pre 7.7147744 0.31215233   0.002      0.002
# 3  Irma vs Pre 0.7382707 0.04690926   0.430      0.430

##time - offshore
seq.sid.lk1o <- data.frame(ps.pm.sid.lk1o@otu_table)
samdf.sid.lk1o <- data.frame(ps.pm.sid.lk1o@sam_data)

dist.sid.lk1o <- vegdist(seq.sid.lk1o)

bet.sid.lk1o <- betadisper(dist.sid.lk1o,samdf.sid.lk1o$time)
#anova(bet.ps.pm)
permutest(bet.sid.lk1o,pairwise=TRUE,permutations=999)#ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.12883 0.064413 0.8942    999  0.433
## Residuals 23 1.65672 0.072031                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.86900 0.310
## Post 0.88312         0.341
## Pre  0.31083 0.33073
plot(bet.sid.lk1o) 

adonis2(seq.sid.lk1o ~ time, data=samdf.sid.lk1o, permutations=999) #p<0.01**
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.lk1o ~ time, data = samdf.sid.lk1o, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)   
## time      2   2.0304 0.27913 4.4531  0.006 **
## Residual 23   5.2434 0.72087                 
## Total    25   7.2738 1.00000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.sid.lk1o, factors=samdf.sid.lk1o$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs  F.Model        R2 p.value p.adjusted
## 1  Post vs Pre 7.621797 0.3226595   0.001      0.001
## 2 Post vs Irma 3.638140 0.1951987   0.021      0.021
## 3  Pre vs Irma 2.032004 0.1193050   0.139      0.139
#post different than other two
#          pairs  F.Model        R2 p.value p.adjusted
# 1  Post vs Pre 7.621797 0.3226595   0.008      0.008
# 2 Post vs Irma 3.638140 0.1951987   0.015      0.015
# 3  Pre vs Irma 2.032004 0.1193050   0.138      0.138

3.2.3 Time - UK1 sids

##time
seq.sid.uk1i <- data.frame(ps.pm.sid.uk1i@otu_table)
samdf.sid.uk1i <- data.frame(ps.pm.sid.uk1i@sam_data)

dist.sid.uk1i <- vegdist(seq.sid.uk1i)

bet.sid.uk1i <- betadisper(dist.sid.uk1i,samdf.sid.uk1i$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk1i,pairwise=TRUE,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.0743 0.037151 0.9067    999  0.434
## Residuals 29 1.1883 0.040975                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.33200 0.684
## Post 0.30524         0.232
## Pre  0.69047 0.20993
plot(bet.sid.uk1i) 

adonis2(seq.sid.uk1i ~ time, data=samdf.sid.uk1i, permutations=999) #p = 0.05 .
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.uk1i ~ time, data = samdf.sid.uk1i, permutations = 999)
##          Df SumOfSqs      R2     F Pr(>F)  
## time      2   0.6774 0.12238 2.022  0.041 *
## Residual 29   4.8576 0.87762               
## Total    31   5.5349 1.00000               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.sid.uk1i, factors=samdf.sid.uk1i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs  F.Model         R2 p.value p.adjusted
## 1  Post vs Pre 2.804997 0.13482324   0.007      0.007
## 2 Post vs Irma 2.497810 0.11102459   0.037      0.037
## 3  Pre vs Irma 0.604099 0.02931936   0.568      0.568
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Post vs Pre 2.804997 0.13482324   0.011      0.011
# 2 Post vs Irma 2.497810 0.11102459   0.030      0.030
# 3  Pre vs Irma 0.604099 0.02931936   0.537      0.537

#post different from other two

##time - offshore
seq.sid.uk1o <- data.frame(ps.pm.sid.uk1o@otu_table)
samdf.sid.uk1o <- data.frame(ps.pm.sid.uk1o@sam_data)

dist.sid.uk1o <- vegdist(seq.sid.uk1o)

bet.sid.uk1o <- betadisper(dist.sid.uk1o,samdf.sid.uk1o$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk1o,pairwise=TRUE,permutations=999)#ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.05342 0.026708 0.3119    999  0.742
## Residuals 30 2.56879 0.085626                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.96000 0.525
## Post 0.96642         0.501
## Pre  0.51116 0.50881
plot(bet.sid.uk1o) 

adonis2(seq.sid.uk1o ~ time, data=samdf.sid.uk1o, permutations=999) #ns
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.uk1o ~ time, data = samdf.sid.uk1o, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)
## time      2   0.7656 0.10339 1.7297  0.125
## Residual 30   6.6393 0.89661              
## Total    32   7.4049 1.00000
pairwise.adonis(seq.sid.uk1o, factors=samdf.sid.uk1o$time, permutations=999) #ns
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs   F.Model         R2 p.value p.adjusted
## 1  Post vs Pre 2.3684810 0.11628167   0.076      0.076
## 2 Post vs Irma 2.3930846 0.10229881   0.081      0.081
## 3  Pre vs Irma 0.5068001 0.02356465   0.465      0.465

3.2.4 Time - UK2 sids

##time - inshore
seq.sid.uk2i <- data.frame(ps.pm.sid.uk2i@otu_table)
samdf.sid.uk2i <- data.frame(ps.pm.sid.uk2i@sam_data)

dist.sid.uk2i <- vegdist(seq.sid.uk2i)

bet.sid.uk2i <- betadisper(dist.sid.uk2i,samdf.sid.uk2i$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk2i,pairwise=TRUE,permutations=999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)   
## Groups     2 0.21882 0.109412 5.6276    999  0.009 **
## Residuals 31 0.60270 0.019442                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          Irma     Post   Pre
## Irma          0.055000 0.049
## Post 0.057523          0.015
## Pre  0.050815 0.013228
##Pre & irma are different
#          Irma     Post   Pre
# Irma          0.439000 0.016
# Post 0.431974          0.187
# Pre  0.017649 0.188058      
plot(bet.sid.uk2i) #pre less constrained

adonis2(seq.sid.uk2i ~ time, data=samdf.sid.uk2i, permutations=999) #p<0.05*
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.uk2i ~ time, data = samdf.sid.uk2i, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## time      2   1.1924 0.17282 3.2385  0.001 ***
## Residual 31   5.7073 0.82718                  
## Total    33   6.8997 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.sid.uk2i, factors=samdf.sid.uk2i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs  F.Model         R2 p.value p.adjusted
## 1 Post vs Irma 3.494325 0.13706285   0.001      0.001
## 2  Post vs Pre 3.981944 0.18114612   0.002      0.002
## 3  Irma vs Pre 1.932622 0.08075262   0.120      0.120
##post different than other two
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Pre vs Post 1.958299 0.05766776   0.028      0.028
# 2  Pre vs Irma 1.398358 0.03461423   0.164      0.164
# 3 Post vs Irma 1.613048 0.03971747   0.050      0.050

##time - offshore
seq.sid.uk2o <- data.frame(ps.pm.sid.uk2o@otu_table)
samdf.sid.uk2o <- data.frame(ps.pm.sid.uk2o@sam_data)

dist.sid.uk2o <- vegdist(seq.sid.uk2o)

bet.sid.uk2o <- betadisper(dist.sid.uk2o,samdf.sid.uk2o$time)
#anova(bet.ps.pm)
permutest(bet.sid.uk2o,pairwise=TRUE,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.08429 0.042143 0.8818    999  0.428
## Residuals 27 1.29039 0.047792                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.31900 0.804
## Post 0.33697         0.110
## Pre  0.77312 0.11717
plot(bet.sid.uk2o) 

adonis2(seq.sid.uk2o ~ time, data=samdf.sid.uk2o, permutations=999) #ns
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.sid.uk2o ~ time, data = samdf.sid.uk2o, permutations = 999)
##          Df SumOfSqs      R2     F Pr(>F)
## time      2   0.6775 0.09085 1.349  0.226
## Residual 27   6.7806 0.90915             
## Total    29   7.4581 1.00000
pairwise.adonis(seq.sid.uk2o, factors=samdf.sid.uk2o$time, permutations=999) #ns
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs   F.Model         R2 p.value p.adjusted
## 1 Post vs Irma 1.6446957 0.08372212   0.116      0.116
## 2  Post vs Pre 2.0931451 0.12245524   0.095      0.095
## 3  Irma vs Pre 0.6640866 0.03065380   0.471      0.471

3.3 Plot by time - sids

3.3.1 PCOA LK1 sids

ps.pm.sid.lk1i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.lk1i)!=0,ps.pm.sid.lk1i)
ord.pm.sid.lk1i.no0 <- ordinate(ps.pm.sid.lk1i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.lk1i.no0, ord.pm.sid.lk1i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

ps.pm.sid.lk1o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.lk1o)!=0,ps.pm.sid.lk1o)
ord.pm.sid.lk1o.no0 <- ordinate(ps.pm.sid.lk1o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.lk1o.no0, ord.pm.sid.lk1o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

3.3.2 PCOA UK1 sids

ps.pm.sid.uk1i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk1i)!=0,ps.pm.sid.uk1i)
ord.pm.sid.uk1i.no0 <- ordinate(ps.pm.sid.uk1i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.uk1i.no0, ord.pm.sid.uk1i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

ps.pm.sid.uk1o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk1o)!=0,ps.pm.sid.uk1o)
ord.pm.sid.uk1o.no0 <- ordinate(ps.pm.sid.uk1o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.uk1o.no0, ord.pm.sid.uk1o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

3.3.3 PCOA UK2 sids

ps.pm.sid.uk2i.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk2i)!=0,ps.pm.sid.uk2i)
ord.pm.sid.uk2i.no0 <- ordinate(ps.pm.sid.uk2i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.uk2i.no0, ord.pm.sid.uk2i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

ps.pm.sid.uk2o.no0 <- prune_taxa(taxa_sums(ps.pm.sid.uk2o)!=0,ps.pm.sid.uk2o)
ord.pm.sid.uk2o.no0 <- ordinate(ps.pm.sid.uk2o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.sid.uk2o.no0, ord.pm.sid.uk2o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

3.4 Stats by time - rads

3.4.1 Setup

ps.pm.rad <- subset_samples(ps.pm,spp=="srad")

##for time stats
ps.pm.rad.lk1 <- subset_samples(ps.pm.rad,site=="LK1")
ps.pm.rad.uk1 <- subset_samples(ps.pm.rad,site=="UK1")
ps.pm.rad.uk2 <- subset_samples(ps.pm.rad,site=="UK2")

ps.pm.rad.lk1i <- subset_samples(ps.pm.rad.lk1,zone=="I")
ps.pm.rad.lk1i <- prune_taxa(taxa_sums(ps.pm.rad.lk1i)!=0,ps.pm.rad.lk1i)
ps.pm.rad.lk1o <- subset_samples(ps.pm.rad.lk1,zone=="O")
ps.pm.rad.lk1o <- prune_taxa(taxa_sums(ps.pm.rad.lk1o)!=0,ps.pm.rad.lk1o)

ps.pm.rad.uk1i <- subset_samples(ps.pm.rad.uk1,zone=="I")
ps.pm.rad.uk1i <- prune_taxa(taxa_sums(ps.pm.rad.uk1i)!=0,ps.pm.rad.uk1i)
ps.pm.rad.uk1o <- subset_samples(ps.pm.rad.uk1,zone=="O")
ps.pm.rad.uk1o <- prune_taxa(taxa_sums(ps.pm.rad.uk1o)!=0,ps.pm.rad.uk1o)

ps.pm.rad.uk2i <- subset_samples(ps.pm.rad.uk2,zone=="I")
ps.pm.rad.uk2i <- prune_taxa(taxa_sums(ps.pm.rad.uk2i)!=0,ps.pm.rad.uk2i)
ps.pm.rad.uk2o <- subset_samples(ps.pm.rad.uk2,zone=="O")
ps.pm.rad.uk2o <- prune_taxa(taxa_sums(ps.pm.rad.uk2o)!=0,ps.pm.rad.uk2o)

3.4.2 Time - LK1 rads

##time - inshore
seq.rad.lk1i <- data.frame(ps.pm.rad.lk1i@otu_table)
samdf.rad.lk1i <- data.frame(ps.pm.rad.lk1i@sam_data)

dist.rad.lk1i <- vegdist(seq.rad.lk1i)

bet.rad.lk1i <- betadisper(dist.rad.lk1i,samdf.rad.lk1i$time)
#anova(bet.ps.pm)
permutest(bet.rad.lk1i,pairwise=TRUE,permutations=999) 
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
## Groups     2 0.90197 0.45098 31.687    999  0.001 ***
## Residuals 21 0.29888 0.01423                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##            Irma       Post   Pre
## Irma            1.0000e-03 0.680
## Post 1.6560e-05            0.001
## Pre  6.7417e-01 4.1092e-05
##Post different than other two
# Response: Distances
#           Df  Sum Sq Mean Sq      F N.Perm Pr(>F)    
# Groups     2 0.90197 0.45098 31.687    999  0.001 ***
# Residuals 21 0.29888 0.01423                         
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Pairwise comparisons:
# (Observed p-value below diagonal, permuted p-value above diagonal)
#            Irma       Post   Pre
# Irma            1.0000e-03 0.664
# Post 1.6560e-05            0.001
# Pre  6.7417e-01 4.1092e-05      
plot(bet.rad.lk1i) 

adonis2(seq.rad.lk1i ~ time, data=samdf.rad.lk1i, permutations=999) #p **
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.rad.lk1i ~ time, data = samdf.rad.lk1i, permutations = 999)
##          Df SumOfSqs     R2      F Pr(>F)   
## time      2   1.1086 0.2633 3.7528  0.007 **
## Residual 21   3.1019 0.7367                 
## Total    23   4.2105 1.0000                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.rad.lk1i, factors=samdf.rad.lk1i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs  F.Model        R2 p.value p.adjusted
## 1  Pre vs Irma 2.630912 0.1798187   0.024      0.024
## 2  Pre vs Post 3.634913 0.1950593   0.033      0.033
## 3 Irma vs Post 4.034110 0.2119411   0.013      0.013
##all different
#          pairs  F.Model        R2 p.value p.adjusted
# 1  Pre vs Irma 2.630912 0.1798187   0.036      0.036
# 2  Pre vs Post 3.634913 0.1950593   0.033      0.033
# 3 Irma vs Post 4.034110 0.2119411   0.024      0.024

##time - offshore
seq.rad.lk1o <- data.frame(ps.pm.rad.lk1o@otu_table)
samdf.rad.lk1o <- data.frame(ps.pm.rad.lk1o@sam_data)

dist.rad.lk1o <- vegdist(seq.rad.lk1o)

bet.rad.lk1o <- betadisper(dist.rad.lk1o,samdf.rad.lk1o$time)
#anova(bet.ps.pm)
permutest(bet.rad.lk1o,pairwise=TRUE,permutations=999)#ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.05765 0.028824 0.3725    999  0.674
## Residuals 18 1.39299 0.077388                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.35000 0.779
## Post 0.34828         0.606
## Pre  0.77610 0.61360
plot(bet.rad.lk1o) 

adonis2(seq.rad.lk1o ~ time, data=samdf.rad.lk1o, permutations=999) #ns
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.rad.lk1o ~ time, data = samdf.rad.lk1o, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)
## time      2   0.5146 0.09471 0.9416  0.453
## Residual 18   4.9186 0.90529              
## Total    20   5.4332 1.00000
pairwise.adonis(seq.rad.lk1o, factors=samdf.rad.lk1o$time, permutations=999)#ns
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs   F.Model         R2 p.value p.adjusted
## 1  Irma vs Pre 0.7328966 0.05755930   0.355      0.355
## 2 Irma vs Post 1.0368385 0.07953144   0.352      0.352
## 3  Pre vs Post 1.0230342 0.07855575   0.303      0.303

3.4.3 Time - UK1 rads

##time
seq.rad.uk1i <- data.frame(ps.pm.rad.uk1i@otu_table)
samdf.rad.uk1i <- data.frame(ps.pm.rad.uk1i@sam_data)

dist.rad.uk1i <- vegdist(seq.rad.uk1i)

bet.rad.uk1i <- betadisper(dist.rad.uk1i,samdf.rad.uk1i$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk1i,pairwise=TRUE,permutations=999) #ns
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.009434 0.0094338 1.1017    999  0.303
## Residuals 10 0.085631 0.0085631                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##        Irma   Pre
## Irma        0.306
## Pre  0.3186
plot(bet.rad.uk1i) 

adonis2(seq.rad.uk1i ~ time, data=samdf.rad.uk1i, permutations=999) #p *
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.rad.uk1i ~ time, data = samdf.rad.uk1i, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)  
## time      1  0.12448 0.25484 3.4199  0.026 *
## Residual 10  0.36399 0.74516                
## Total    11  0.48847 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.rad.uk1i, factors=samdf.rad.uk1i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
##         pairs F.Model        R2 p.value p.adjusted
## 1 Pre vs Irma 3.41993 0.2548396   0.033      0.033
# adonis2(formula = seq.rad.uk1i ~ time, data = samdf.rad.uk1i, permutations = 999)
#          Df SumOfSqs      R2      F Pr(>F)  
# time      1  0.12448 0.25484 3.4199   0.03 *
# Residual 10  0.36399 0.74516                
# Total    11  0.48847 1.00000

#offshore not available

3.4.4 Time - UK2 rads

##time - inshore
seq.rad.uk2i <- data.frame(ps.pm.rad.uk2i@otu_table)
samdf.rad.uk2i <- data.frame(ps.pm.rad.uk2i@sam_data)

dist.rad.uk2i <- vegdist(seq.rad.uk2i)

bet.rad.uk2i <- betadisper(dist.rad.uk2i,samdf.rad.uk2i$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk2i,pairwise=TRUE,permutations=999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq     F N.Perm Pr(>F)  
## Groups     2 0.22544 0.112720 5.305    999  0.015 *
## Residuals 21 0.44621 0.021248                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##          Irma     Post   Pre
## Irma          0.039000 0.225
## Post 0.033448          0.018
## Pre  0.222245 0.023323
##Pre & irma are different
#          Irma     Post   Pre
# Irma          0.026000 0.218
# Post 0.033448          0.016
# Pre  0.222245 0.023323        
plot(bet.rad.uk2i) #post less constrained

adonis2(seq.rad.uk2i ~ time, data=samdf.rad.uk2i, permutations=999) #p **
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.rad.uk2i ~ time, data = samdf.rad.uk2i, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)    
## time      2   0.8178 0.21125 2.8123  0.001 ***
## Residual 21   3.0534 0.78875                  
## Total    23   3.8712 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pairwise.adonis(seq.rad.uk2i, factors=samdf.rad.uk2i$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs  F.Model         R2 p.value p.adjusted
## 1  Pre vs Irma 1.514262 0.09169421   0.150      0.150
## 2  Pre vs Post 2.938851 0.19672535   0.008      0.008
## 3 Irma vs Post 3.344836 0.18233122   0.002      0.002
##post different than other two
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Pre vs Irma 1.514262 0.09169421   0.161      0.161
# 2  Pre vs Post 2.938851 0.19672535   0.004      0.004
# 3 Irma vs Post 3.344836 0.18233122   0.005      0.005

##time - offshore
seq.rad.uk2o <- data.frame(ps.pm.rad.uk2o@otu_table)
samdf.rad.uk2o <- data.frame(ps.pm.rad.uk2o@sam_data)

dist.rad.uk2o <- vegdist(seq.rad.uk2o)

bet.rad.uk2o <- betadisper(dist.rad.uk2o,samdf.rad.uk2o$time)
#anova(bet.ps.pm)
permutest(bet.rad.uk2o,pairwise=TRUE,permutations=999)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
## Groups     2 0.03579 0.017893 0.2306    999  0.805
## Residuals 23 1.78482 0.077601                     
## 
## Pairwise comparisons:
## (Observed p-value below diagonal, permuted p-value above diagonal)
##         Irma    Post   Pre
## Irma         0.65000 0.562
## Post 0.66724         0.750
## Pre  0.56121 0.74220
##Pre & irma are different
#          Irma     Post   Pre
# Irma          0.026000 0.218
# Post 0.033448          0.016
# Pre  0.222245 0.023323        
plot(bet.rad.uk2o) #post less constrained

adonis2(seq.rad.uk2o ~ time, data=samdf.rad.uk2o, permutations=999) #p **
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = seq.rad.uk2o ~ time, data = samdf.rad.uk2o, permutations = 999)
##          Df SumOfSqs      R2      F Pr(>F)
## time      2   0.7432 0.10812 1.3941  0.178
## Residual 23   6.1308 0.89188              
## Total    25   6.8740 1.00000
pairwise.adonis(seq.rad.uk2o, factors=samdf.rad.uk2o$time, permutations=999)
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
## 'adonis' will be deprecated: use 'adonis2' instead
##          pairs   F.Model         R2 p.value p.adjusted
## 1  Post vs Pre 2.0874958 0.10936457   0.081      0.081
## 2 Post vs Irma 1.6509889 0.10548783   0.118      0.118
## 3  Pre vs Irma 0.4521508 0.02926135   0.682      0.682
##post different than other two
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Pre vs Irma 1.514262 0.09169421   0.161      0.161
# 2  Pre vs Post 2.938851 0.19672535   0.004      0.004
# 3 Irma vs Post 3.344836 0.18233122   0.005      0.005

3.5 Plot by time - rads

3.5.1 PCOA LK1 rads

ps.pm.rad.lk1i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.lk1i)!=0,ps.pm.rad.lk1i)
ord.pm.rad.lk1i.no0 <- ordinate(ps.pm.rad.lk1i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.lk1i.no0, ord.pm.rad.lk1i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

ps.pm.rad.lk1o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.lk1o)!=0,ps.pm.rad.lk1o)
ord.pm.rad.lk1o.no0 <- ordinate(ps.pm.rad.lk1o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.lk1o.no0, ord.pm.rad.lk1o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

3.5.2 PCOA UK1 rads

ps.pm.rad.uk1i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk1i)!=0,ps.pm.rad.uk1i)
ord.pm.rad.uk1i.no0 <- ordinate(ps.pm.rad.uk1i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.uk1i.no0, ord.pm.rad.uk1i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

ps.pm.rad.uk1o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk1o)!=0,ps.pm.rad.uk1o)
ord.pm.rad.uk1o.no0 <- ordinate(ps.pm.rad.uk1o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.uk1o.no0, ord.pm.rad.uk1o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

3.5.3 PCOA UK2 rads

ps.pm.rad.uk2i.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk2i)!=0,ps.pm.rad.uk2i)
ord.pm.rad.uk2i.no0 <- ordinate(ps.pm.rad.uk2i.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.uk2i.no0, ord.pm.rad.uk2i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()

ps.pm.rad.uk2o.no0 <- prune_taxa(taxa_sums(ps.pm.rad.uk2o)!=0,ps.pm.rad.uk2o)
ord.pm.rad.uk2o.no0 <- ordinate(ps.pm.rad.uk2o.no0,"PCoA",distance="bray")

plot_ordination(ps.pm.rad.uk2o.no0, ord.pm.rad.uk2o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  #facet_grid(site~zone)+
  theme_cowplot()
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

4 Script graveyard

Type profile business

setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
#ps <- readRDS("ps.its2.rds")
#ps.rel <- transform_sample_counts(ps, function(x) x / sum(x))
#saveRDS(ps.rel,"ps.its2.rel.rds")
ps1 <- readRDS("ps.its2.rel.rds")
ps <- subset_samples(ps1,site!="UK3")

ps.sid <- subset_samples(ps,spp=="ssid")

ps.sid.pre <- subset_samples(ps.sid,time=="Pre")
ps.sid.irm <- subset_samples(ps.sid,time=="Irma")
ps.sid.pos <- subset_samples(ps.sid,time=="Post")

ps.sid.pre.lk1 <- subset_samples(ps.sid.pre,site=="LK1")
ps.sid.pre.uk1 <- subset_samples(ps.sid.pre,site=="UK1")
ps.sid.pre.uk2 <- subset_samples(ps.sid.pre,site=="UK2")

ps.sid.irm.lk1 <- subset_samples(ps.sid.irm,site=="LK1")
ps.sid.irm.uk1 <- subset_samples(ps.sid.irm,site=="UK1")
ps.sid.irm.uk2 <- subset_samples(ps.sid.irm,site=="UK2")

ps.sid.pos.lk1 <- subset_samples(ps.sid.pos,site=="LK1")
ps.sid.pos.uk1 <- subset_samples(ps.sid.pos,site=="UK1")
ps.sid.pos.uk2 <- subset_samples(ps.sid.pos,site=="UK2")

##for time stats
ps.sid.lk1 <- subset_samples(ps.sid,site=="LK1")
ps.sid.uk1 <- subset_samples(ps.sid,site=="UK1")
ps.sid.uk2 <- subset_samples(ps.sid,site=="UK2")

ps.sid.lk1i <- subset_samples(ps.sid.lk1,zone=="I")
ps.sid.lk1o <- subset_samples(ps.sid.lk1,zone=="O")

ps.sid.uk1i <- subset_samples(ps.sid.uk1,zone=="I")
ps.sid.uk1o <- subset_samples(ps.sid.uk1,zone=="O")

ps.sid.uk2i <- subset_samples(ps.sid.uk2,zone=="I")
ps.sid.uk2o <- subset_samples(ps.sid.uk2,zone=="O")

4.0.1 Plots by time

ps.sid.lk1i.no0 <- prune_taxa(taxa_sums(ps.sid.lk1i)!=0,ps.sid.lk1i)

ord.sid.lk1i.no0 <- ordinate(ps.sid.lk1i.no0,"PCoA",distance="bray")

plot_ordination(ps.sid.lk1i.no0, ord.sid.lk1i.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  theme_cowplot()

ps.sid.lk1o.no0 <- prune_taxa(taxa_sums(ps.sid.lk1o)!=0,ps.sid.lk1o)

ord.sid.lk1o.no0 <- ordinate(ps.sid.lk1o.no0,"PCoA",distance="bray")

plot_ordination(ps.sid.lk1o.no0, ord.sid.lk1o.no0, color ="time")+
  geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  theme_cowplot()

4.0.2 Reef zones

ord.sid.pre.lk1 <- ordinate(ps.sid.pre.lk1,"PCoA",distance="bray")
plot_ordination(ps.sid.pre.lk1, ord.sid.pre.lk1, color ="zone")+
  #geom_point(alpha=0.5)+
  #scale_color_manual(values=c("darkslategray3","darkslategray4","#000004"))+
  #scale_shape_manual(values=c(16,15))+
  stat_ellipse()+
  theme_cowplot()
  #facet_grid(site~zone)

#ggsave(filename="pcoa.types.site.pdf")

5 Stats - type profiles

5.1 Setup

library(vegan)
#remotes::install_github("Jtrachsel/funfuns")
library(funfuns)
library(dplyr)
library(edgeR)

setwd("/Volumes/GoogleDrive/My Drive/Flirma/flirma_networks/fl_its2")
ps <- readRDS("ps.its2.rds")

ps <- subset_samples(ps,site!="UK3")

5.2 Subset

ps.sid <- subset_samples(ps,spp="ssid")
ps.sid <- transform_sample_counts(ps.sid, function(x) x / sum(x))

##for in v off stats
ps.sid.pre <- subset_samples(ps.sid,time=="Pre")
ps.sid.irm <- subset_samples(ps.sid,time=="Irma")
ps.sid.pos <- subset_samples(ps.sid,time=="Post")

ps.sid.pre.lk1 <- subset_samples(ps.sid.pre,site=="LK1")
ps.sid.pre.uk1 <- subset_samples(ps.sid.pre,site=="UK1")
ps.sid.pre.uk2 <- subset_samples(ps.sid.pre,site=="UK2")

ps.sid.irm.lk1 <- subset_samples(ps.sid.irm,site=="LK1")
ps.sid.irm.uk1 <- subset_samples(ps.sid.irm,site=="UK1")
ps.sid.irm.uk2 <- subset_samples(ps.sid.irm,site=="UK2")

ps.sid.pos.lk1 <- subset_samples(ps.sid.pos,site=="LK1")
ps.sid.pos.uk1 <- subset_samples(ps.sid.pos,site=="UK1")
ps.sid.pos.uk2 <- subset_samples(ps.sid.pos,site=="UK2")

##for time stats
ps.sid.lk1 <- subset_samples(ps.sid,site=="LK1")
ps.sid.uk1 <- subset_samples(ps.sid,site=="UK1")
ps.sid.uk2 <- subset_samples(ps.sid,site=="UK2")

ps.sid.lk1i <- subset_samples(ps.sid.lk1,zone=="I")
ps.sid.lk1o <- subset_samples(ps.sid.lk1,zone=="O")

ps.sid.uk1i <- subset_samples(ps.sid.uk1,zone=="I")
ps.sid.uk1o <- subset_samples(ps.sid.uk1,zone=="O")

ps.sid.uk2i <- subset_samples(ps.sid.uk2,zone=="I")
ps.sid.uk2o <- subset_samples(ps.sid.uk2,zone=="O")

5.3 Adonisii - sids

5.3.1 Reef zones - LK1

##reef zone pre
seq.sid.pre.lk1 <- data.frame(ps.sid.pre.lk1@otu_table)
samdf.sid.pre.lk1 <- data.frame(ps.sid.pre.lk1@sam_data)

dist.sid.pre.lk1 <- vegdist(seq.sid.pre.lk1)
##zone
bet.sid.pre.lk1 <- betadisper(dist.sid.pre.lk1,samdf.sid.pre.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.pre.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pre.lk1) #ns

adonis2(seq.sid.pre.lk1 ~ zone, data=samdf.sid.pre.lk1, permutations=999) #ns

##reef zone irm
seq.sid.irm.lk1 <- data.frame(ps.sid.irm.lk1@otu_table)
samdf.sid.irm.lk1 <- data.frame(ps.sid.irm.lk1@sam_data)

dist.sid.irm.lk1 <- vegdist(seq.sid.irm.lk1)
##zone
bet.sid.irm.lk1 <- betadisper(dist.sid.irm.lk1,samdf.sid.irm.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.irm.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.irm.lk1) #ns

adonis2(seq.sid.irm.lk1 ~ zone, data=samdf.sid.irm.lk1, permutations=999)
#adonis2(formula = seq.sid.irm.lk1 ~ zone, data = samdf.sid.irm.lk1, permutations = 999)
#          Df SumOfSqs      R2     F Pr(>F)  
# zone      1   0.8368 0.06787 1.966  0.038 *
# Residual 27  11.4929 0.93213               
# Total    28  12.3298 1.00000               
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

##reef zone post
seq.sid.pos.lk1 <- data.frame(ps.sid.pos.lk1@otu_table)
samdf.sid.pos.lk1 <- data.frame(ps.sid.pos.lk1@sam_data)

dist.sid.pos.lk1 <- vegdist(seq.sid.pos.lk1)
##zone
bet.sid.pos.lk1 <- betadisper(dist.sid.pos.lk1,samdf.sid.pos.lk1$zone)
#anova(bet.ps)
permutest(bet.sid.pos.lk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pos.lk1) #ns

adonis2(seq.sid.pos.lk1 ~ zone, data=samdf.sid.pos.lk1, permutations=999) #ns

5.3.2 Reef zones - UK1

##reef zone pre
seq.sid.pre.uk1 <- data.frame(ps.sid.pre.uk1@otu_table)
samdf.sid.pre.uk1 <- data.frame(ps.sid.pre.uk1@sam_data)

dist.sid.pre.uk1 <- vegdist(seq.sid.pre.uk1)
##zone
bet.sid.pre.uk1 <- betadisper(dist.sid.pre.uk1,samdf.sid.pre.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.pre.uk1,pairwise=FALSE,permutations=999)
plot(bet.sid.pre.uk1)
#          Df SumOfSqs      R2      F Pr(>F)  
# zone      1    0.914 0.06604 2.0507  0.023 *
# Residual 29   12.925 0.93396                
# Total    30   13.839 1.00000                

adonis2(seq.sid.pre.uk1 ~ zone, data=samdf.sid.pre.uk1, permutations=999) 
# adonis2(formula = seq.sid.pre.uk1 ~ zone, data = samdf.sid.pre.uk1, permutations = 999)
#          Df SumOfSqs      R2      F Pr(>F)  
# zone      1    0.914 0.06604 2.0507  0.021 *
# Residual 29   12.925 0.93396                
# Total    30   13.839 1.00000                
# ---

##reef zone irm
seq.sid.irm.uk1 <- data.frame(ps.sid.irm.uk1@otu_table)
samdf.sid.irm.uk1 <- data.frame(ps.sid.irm.uk1@sam_data)

dist.sid.irm.uk1 <- vegdist(seq.sid.irm.uk1)
##zone
bet.sid.irm.uk1 <- betadisper(dist.sid.irm.uk1,samdf.sid.irm.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.irm.uk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.irm.uk1) #ns

adonis2(seq.sid.irm.uk1 ~ zone, data=samdf.sid.irm.uk1, permutations=999)
# adonis2(formula = seq.sid.irm.uk1 ~ zone, data = samdf.sid.irm.uk1, permutations = 999)
#          Df SumOfSqs      R2      F Pr(>F)   
# zone      1   1.3958 0.09887 3.2915  0.002 **
# Residual 30  12.7216 0.90113                 
# Total    31  14.1174 1.00000          

##reef zone post
seq.sid.pos.uk1 <- data.frame(ps.sid.pos.uk1@otu_table)
samdf.sid.pos.uk1 <- data.frame(ps.sid.pos.uk1@sam_data)

dist.sid.pos.uk1 <- vegdist(seq.sid.pos.uk1)
##zone
bet.sid.pos.uk1 <- betadisper(dist.sid.pos.uk1,samdf.sid.pos.uk1$zone)
#anova(bet.ps)
permutest(bet.sid.pos.uk1,pairwise=FALSE,permutations=999) #ns
plot(bet.sid.pos.uk1) #ns

adonis2(seq.sid.pos.uk1 ~ zone, data=samdf.sid.pos.uk1, permutations=999) #ns

5.3.3 Reef zone - UK2

5.3.4 Time - LK1 sids

##time - inshore
seq.sid.lk1i <- data.frame(ps.sid.lk1i@otu_table)
samdf.sid.lk1i <- data.frame(ps.sid.lk1i@sam_data)

dist.sid.lk1i <- vegdist(seq.sid.lk1i)

bet.sid.lk1i <- betadisper(dist.sid.lk1i,samdf.sid.lk1i$time)
#anova(bet.ps)
permutest(bet.sid.lk1i,pairwise=TRUE,permutations=999)
#post different than other two, looks more constrained
# Permutation test for homogeneity of multivariate dispersions
# Permutation: free
# Number of permutations: 999
# 
# Response: Distances
#           Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
# Groups     2 0.05045 0.0252251 4.4942    999  0.016 *
# Residuals 46 0.25819 0.0056128                       
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Pairwise comparisons:
# (Observed p-value below diagonal, permuted p-value above diagonal)
#           Irma      Post   Pre
# Irma           0.0240000 0.349
# Post 0.0239826           0.006
# Pre  0.3359886 0.0082901  
plot(bet.sid.lk1i) 

adonis2(seq.sid.lk1i ~ time, data=samdf.sid.lk1i, permutations=999) #p<0.001***

pairwise.adonis(seq.sid.lk1i, factors=samdf.sid.lk1i$time, permutations=999)
##pre diff than post
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Pre vs Irma 1.273944 0.04208053   0.272      0.272
# 2  Pre vs Post 3.211953 0.08869871   0.001      0.001
# 3 Irma vs Post 1.507718 0.04785234   0.124      0.124

##time - offshore
seq.sid.lk1o <- data.frame(ps.sid.lk1o@otu_table)
samdf.sid.lk1o <- data.frame(ps.sid.lk1o@sam_data)

dist.sid.lk1o <- vegdist(seq.sid.lk1o)

bet.sid.lk1o <- betadisper(dist.sid.lk1o,samdf.sid.lk1o$time)
#anova(bet.ps)
permutest(bet.sid.lk1o,pairwise=TRUE,permutations=999)#ns
plot(bet.sid.lk1o) 

adonis2(seq.sid.lk1o ~ time, data=samdf.sid.lk1o, permutations=999) #p<0.01**

pairwise.adonis(seq.sid.lk1o, factors=samdf.sid.lk1o$time, permutations=999)
#post different than other two
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Irma vs Pre 1.549658 0.05072586   0.132      0.132
# 2 Irma vs Post 2.674018 0.08442309   0.001      0.001
# 3  Pre vs Post 2.582747 0.07926733   0.014      0.014

5.3.5 Time - UK1 sids

##time
seq.sid.uk1i <- data.frame(ps.sid.uk1i@otu_table)
samdf.sid.uk1i <- data.frame(ps.sid.uk1i@sam_data)

dist.sid.uk1i <- vegdist(seq.sid.uk1i)

bet.sid.uk1i <- betadisper(dist.sid.uk1i,samdf.sid.uk1i$time)
#anova(bet.ps)
permutest(bet.sid.uk1i,pairwise=TRUE,permutations=999) #ns
plot(bet.sid.uk1i) 

adonis2(seq.sid.uk1i ~ time, data=samdf.sid.uk1i, permutations=999) #p<0.01**

pairwise.adonis(seq.sid.uk1i, factors=samdf.sid.uk1i$time, permutations=999)
#         pairs  F.Model         R2 p.value p.adjusted
# 1  Post vs Pre 1.903144 0.07347155   0.026      0.026
# 2 Post vs Irma 2.306884 0.08149552   0.022      0.022
# 3  Pre vs Irma 1.898341 0.05600101   0.047      0.047

#all different

5.3.6 Time - UK2 sids

##time
seq.sid.uk2i <- data.frame(ps.sid.uk2i@otu_table)
samdf.sid.uk2i <- data.frame(ps.sid.uk2i@sam_data)

dist.sid.uk2i <- vegdist(seq.sid.uk2i)

bet.sid.uk2i <- betadisper(dist.sid.uk2i,samdf.sid.uk2i$time)
#anova(bet.ps)
permutest(bet.sid.uk2i,pairwise=TRUE,permutations=999)
##Pre & irma are different
#          Irma     Post   Pre
# Irma          0.439000 0.016
# Post 0.431974          0.187
# Pre  0.017649 0.188058      
plot(bet.sid.uk2i) #pre less constrained

adonis2(seq.sid.uk2i ~ time, data=samdf.sid.uk2i, permutations=999) #p<0.05*

pairwise.adonis(seq.sid.uk2i, factors=samdf.sid.uk2i$time, permutations=999)
##post different than other two
#          pairs  F.Model         R2 p.value p.adjusted
# 1  Pre vs Post 1.958299 0.05766776   0.028      0.028
# 2  Pre vs Irma 1.398358 0.03461423   0.164      0.164
# 3 Post vs Irma 1.613048 0.03971747   0.050      0.050
ps.sid.rel <- transform_sample_counts(ps.sid, function(x) x / sum(x))
## Presence/absence

sppAbun <- as.data.frame(ps.sid.rel@otu_table)

sppAbun[sppAbun > 0] <- 1 #converts from abundance to P/A
sppAbun

rowSums(sppAbun)

morethan1 <- sppAbun[rowSums(sppAbun)>1,]

find.top.taxa <- function(x,taxa){
  require(phyloseq)
  top.taxa <- tax_glom(x, taxa)
  otu <- otu_table(t(top.taxa)) # remove the transformation if using a merge_sample object
  tax <- tax_table(top.taxa)
  j<-apply(otu,1,which.max)
  k <- j[!duplicated(j)]
  l <- data.frame(tax[k,])
  m <- data.frame(otu[,k])
  s <- as.name(taxa)
  colnames(m) = l[,taxa]
  n <- colnames(m)[apply(m,1,which.max)]
  m[,taxa] <- n
  return(m)
}
find.top.taxa(top.pdy,"Phylum")

#chi-square, top vs side
chisq.test(SymPortalDataK_1$variable, SymPortalDataK_1$label, correct=FALSE)
#Pearson's Chi-squared test

#data:  SymPortalDataK_1$variable and SymPortalDataK_1$label
#X-squared = 5.5824, df = 6, p-value = 0.4716

5.4 Find “most frequent” sym

Just ran once to get the info

find.top.asv <- function(x,num){
  require(phyloseq)
  require(magrittr)
  
  otu <- otu_table(x)
  tax <- tax_table(x)
  j1 <- apply(otu,1,sort,index.return=T, decreasing=T) # modifying which.max to return a list of sorted index
  j2 <- lapply(j1,'[[',"ix") # select for index
  
  l <- data.frame(unique(tax@.Data[unlist(j2),]))
  m <- data.frame(otu@.Data[,unique(unlist(j2))])
  n <- apply(m,1,sort,index.return=T, decreasing=T) %>%
    lapply('[[',"ix") %>%  # Extract index
    lapply(head,n=num) # This to returns the top x tax

  p <- list()
  for(i in 1:length(n)){
    p[[i]]<- colnames(m)[n[[i]]]
  }
  m$taxa <- p
  return(m)
}

top.asvs <- find.top.asv(ps.sid,5)
top.asvs$taxa <- as.character(top.asvs$taxa)

#write.csv(top.asvs,file="top_sym.csv",row.names=TRUE)
table <- read.csv("silly_its_table.csv",header=T)

heatmap(as.matrix(table[,2:11]), labRow=table[,1],rowv=NA)

6 Flirma scripts scraps

6.1 Collapsing type profiles via phylogeny tree

Started with file: “188_20211214_03_DBV_20211215T030311.seqs.fasta” in the post-med seqs folder. Also started with the taxa info ‘Majority_ITS2_sequence", removed the slashes between types, & removed duplicates. That left 24 unique ITS2 types, put these in ’names.txt’. Want to extract the relevant sequences from the fasta file now

for word in $(cat names.txt); do echo $word; grep -w -A 1 $word postmedseqs.fasta >> output.fasta; done

Important notes: produced a fasta file with ‘–’ in between some desired entries, removed that manually, but I’m sure there’s a better way than that. It also gave me an unwanted one (A4.3), took that out manually. Easy to do with 24 types, but would be harder with more than that.. processed file was named output_types_cleaned.fasta

Then went to ngphylogeny.fr to plot.

More notes: - A is fine, the majority sequences are pretty disparate. - The Bs are a total mess, will plot below - Cs are fine - Ds are fine

#facet_grid(site~zone)
# 
# #devtools::install_github("david-barnett/microViz@0.9.3") # check 0.9.3 is the latest release
# library("microViz")
# 
# ps.sid.lk1i.no0 %>%
#   ps_mutate(
#     time = as.numeric(time == "time")
#   ) %>%
#   ord_calc(
#     constraints = c("time"),
#     # method = "RDA", # Note: you can specify RDA explicitly, and it is good practice to do so, but microViz can guess automatically that you want an RDA here (helpful if you don't remember the name?)
#     scale_cc = FALSE # doesn't make a difference
#   ) %>%
#   plot_ord(
#     colour = "time", size = 2, alpha = 0.5, shape = "active",
#     plot_taxa = 1:8
#   )